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problem  for  an  autonomous  system  of  N ordinary  differential  equations.  Using  fast 
power  series  techniques,  we  exhibit  an  algorithm  for  the  p^^-order  Taylor  series 
method  requiring  only  0(p^  In  p)  arithmetic  operations  per  step  as  p -»  ♦».  (Moreover, 
we  show  that  any  such  algorithm  requires  at  least  CKp^)  operations  per  step.)  We 
compute  the  order  which  minimizes  the  complexity  bounds  for  Taylor  series  and  linear 
Runge-Kutta  methods,  and  show  that  in  all  cases,  this  optimal  order  increases  as  the 
error  criterion  i decreases,  tending  to  infinity  as  t ten^  to  rero.  Finally,  we  show 
that  if  certain  derivatives  are  easy  to  evaluate,  then  Taylor  series  methods  are 
asymptotically  better  than  linear  Runge*Kutta  methods  for  problems  of  small 
dimension  N. 
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1.  Introduction 


Let  be  a set  of  points  in  the  real  N-ditnensional  linear  space  R 


such  that  the  initial  value  orobletn  of  finding  a function 


be  a set  of  operators  on  R 


X : [0,  1]  -»  R*'  satisfying 


has  a unique  solution  for  every  (xq  , v)  f ; we  assume  that  x is  analytic  on  [0,  1]^ 
The  autonomous  form  of  this  system  is  no  restriction,  since  any  non-autonomous 


system  may  be  made  autonomous  by  increasing  the  dimension  of  the  system  by  one. 

In  Werschulz  [76],  we  looked  at  the  computational  complexity  of  using  one-step 
methods  to  generate  an  approximate  solution  to  (1.1)  on  an  equidistant  grid  in  the 
sense  of  Stetter  [73];  that  is,  the  methods  considered  computed  approximations  X{  to 
x(ih)  by  the  recursion 

(1.2)  Xj+j  - Xj  + h <?(Xj , h)  (0  S i S n - 1 , n - h"^) , 

where  h ■ n’^  is  the  step-size  of  a grid  with  n points,  and  p is  the  increment  function 
(Henrici  [62])  for  the  method.  (To  be  brief,  we  will  refer  to  "the  method  p.")  In  that 
paper,  we  discussed  the  problem  of  optimal  order  and  minimal  complexity  for  rather 


general  classes  of  one-step  methods. 

In  this  paper,  we  will  use  the  techniques  and  results  of  Werschulz  [76]  to 
analyze  the  complexity  of  using  Taylor  series  methods  and  linear  Runge-Kutta  methods 
to  generate  approximate  solutions  whose  error  does  not  exceed  i.  The  model  of  v 
computation,  error  measure,  and  complexity  measure  to  be  used  are  described  in 
Section  2,  as  well  as  the  relevant  results  from  Werschulz  [76].  4 
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We  discuss  the  complexity  of  Taylor  series  methods  in  Section  3.  Using  the  feat 
power  series  techniques  of  Brent  and  Kung  [761  show  that  (Xp^  In  p)  erithmetle 
operations  suffice  to  compute  the  p^^-order  Taylor  series  approximation)  moreoveri  we 
show  that  CKp^)  operations  are  necessary.  In  Section  4,  we  discuss  the  complexity  of 
linear  Runge-Kutta  methods.  In  both  Sections,  we  compute  lower  and  upper  bounds  on 
the  compiexity  using  a fixed  method  of  given  order)  these  results  are  then  used  to 
compute  optimal  orders  which  minimize  these  compiexity  bounds.  We  show  that  in  all 
cases,  the  optimal  order  increases  as  i decreases,  tending  to  infinity  as  a tends  to 
zero. 

Finally,  we  compare  these  two  classes  of  methods  in  Section  5,  where  we  show 
that  if  the  partial  derivatives  of  v are  easy  to  evaluate,  then  Taylor  series  methods  are 
asymptotically  better  (as  • tends  to  zero)  than  linear  Runge-Kutta  methods  for 
problems  of  small  dimension  N. 
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2.  Preliminary  Results 

Before  proceeding  any  further,  we  will  establish  some  notational  conventions. 
Let  X be  an  ordered  ring;  then  X'*'  and  X'*"*^  respectively  denote  the  nonnegative  and 
positive  elements  of  X (This  is  used  in  the  cases  X ■ It,  the  real  numbers,  and 
X ■>  Z,  the  integers.)  The  symbol  means  "is  defined  to  be,"  while  "ii"  means  "Is 
identically  equal  to."  We  use  "I"  to  denote  the  unit  interval  [0, 1].  The  symbol  "V"  is 
used  to  denote  the.  gradient  of  a mapping.  The  notations  "x  I a"  and  "x  t a"  are  used 
to  indicate  one-sided  limits,  as  in  Buck  [65].  Finally,  we  write  "(a.b)^”  to  indicate  the 
c^^  part  of  equation  (a.b),  as  in  Gurtin  [75]. 

We  next  describe  the  model  of  computation  to  be  used.  We  assume  only  that  ail 
arithmetic  operations  are  performed  exactly  in  R (i.e.,  infinite -precision  arithmetic)  and 
that  for  any  algorithm  to  be  considered  for  the  solution  of  (1.1),  a set  of  procedures  is 
given  for  the  computation  of  any  information  about  v required  by  that  algorithm.  (For 
instance,  with  Runge-Kutta  methods,  we  must  be  able  to  compute  v at  any  point  in  its 
domain.) 

In  addition,  we  must  pick  an  error  measure,  so  that  we  may  measure  the 
discrepancy  between  the  approximate  solution  produced  by  f (via  (1.2))  and  the  true 
solution.  For  the  sake  of  definiteness,  we  use  the  global  error 

(2.1)  aG^e.h)  :•  max  o^j^n  IW'W  - Xj||  , 

where  ||  ’ ||  is  a norm  on  R^ . Other  error  measures  may  be  used,  such  as  the  local 
error  per  step  a|^  and  the  local  error  per  unit  a|_y  (see  Hanrici  [62]  and 
Stetter  [73]  for  definitions)}  this  would  involve  only  a slight  modification  of  the  results 
contained  in  the  sequel. 


Wa  finally  describe  the  complexity  measure  to  be  used.  Let  4 « * P ^ Z 


be  a basic  sequence  in  the  sense  of  Werschulz  [76];  that  is,  there  exist  functions 
K ! -♦  !?■*■  and  xi  , xii : R*  -*  R*  such  that 


where 


(2.3)  0 < »l(p>  a «(p,h)  S mjip)  < *00  for  h < I . 

We  say  that  has  order  p.  This  is  a slight  extension  of  the  definition  of  order  given 


in  Cooper  and  Verner  [72]i  the  function  introduced  here  is  necessary  and  sufficient 
for  the  "order"  of  a method  to  be  unique.  (For  the  sake  of  exposition,  we  assume  that 
*L  and  «y  are  analytic  on  R'*',  and  that  lim  p^o  «l(p)^^P  and  lim  «u(p)^^P  exist  and 
are  positive  real  numbers;  this  will  always  be  the  case  in  the  examples  we  consider.) 
Then  we  will  be  interested  jn  the  total  number  gf  arithmetic  operations  C(p,«)  required 


(2.4)  ^ ‘ ® * ' 

for  a given  p and  a given  a.  (Here  e is  the  base  of  the  natural  logarithms.)  We 


suppose  that  0 < t < 1,  so  that  a is  positive.  Clearly,  a increases  as  • decreases,  and  « 


tends  to  infinity  as  « tends  to  zero. 


In  the  methods  we  consider,  we  may  write 


where  n is  the  minimal  number  of  steps  required  and  the  cost  per  step  c(p)  is  the 
number  of  arithmetic  operations  required  tor  the  method  of  order  p.  As  in  Traub  and 
WoiniaKowsKi  [76],  we  shall  express  the  cost  per  step  associated  with  pp  in  the  form 
(2.6)  c(p)  e(9tp(v))  ♦ d(p)  . 

Here  Wp(v)  is  the  information  about  v required  to  perform  one  step  of  ^p,  and  we 
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write  e(Kp(v»  for  the  informationel  cQsi  of  we  call  d(p)  the  combinatory  coat 


■X 

I • 


of  Wp. 

Note  that  we  explicitly  indicate  the  dependence  of  f|p  on  v,  to  that  we  may 
compare  the  cost  of  <say)  an  evaluation  of  v with  a scalar  arithmetic  operation. 
Basically,  e(92p(v»  measures  the  cost  of  getting  new  data  about  v required  by  ^p, 
while  d(p)  measures  the  cost  of  combining  this  new  data  to  get  an  approximate  value 
of  the  solution  at  a new  point.  For  example,  Euler’s  method  in 

*i+l  “ ><i  ♦ v(Xj) 

has  informational  cost  z|^,  e(vj) , where  , ... , Vf^  are  the  components  of  and  for 
any  function  w:  •*  R,  we  define 

(2.7)  p(«)  :•<  cost  of  evaluating  w at  one  point  . 

The  combinatory  cost  is  2rj  arithmetic  operations,  i.e.,  one  scalar  multiplication  and  one 
scalar  addition  for  each  of  the  N components. 

We  must  now  face  a problem  that  occurs  in  almost  alt  areas  of  complexity 
theory.  The  number  of  operations  c(p)  required  for  one  step  of  a p*^-order  method  is 
usually  unknown  psr  sg;  we  only  have  bounds  of  the  form 

(2.8)  CL<p)  S c(p)  i cj(p)  . 

That  is,  C|_(p)  is  a lower  bound  on  the  number  of  operations  required  per  step,  usually 
derived  via  theoretical  considerations,  and  Cy(p)  is  an  upper  bound  on  the  number  of 
operations  required  per  step,  which  is  derived  by  exhibiting  an  algorithm  ‘for 
computing  the  p^^-order  method.  (In  what  follows,  we  shall  assume  that  the  functions 
C|_ , Cy  : R'*'  ->  R'*’  are  analytic,  although  this  requirement  may  be  greatly  weakened. 
However,  this  assumption  holds  for  all  examples  that  we  consider.) 

From  the  discussion  in  Section  3 of  Werschulz  [76],  we  find  that  the  step-size  h 
must  satisfy 


■'  1 
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• i 

(2.9)  hu(p,«)  S h S hL(p^) , 

i 

wherp 

(2.10)  h|^(p#)  *l(p)“^^P  e"*/p  and  hu(p,«)  «ij(p)"^^P  • i 

Using  (2.5),  (2.8),  (2.9),  and  (2.10),  we  may  find  bounds  on  the  complexity  C(p,e).  j 

i 

Theorem  2.1;  Define  \ -i 

CL(p,e)  fL(p)e*/P,  where  f^fp)  «l^p)^^P  Cl(p)  , 

and 

Cu(p,«)  f|j(p)  e“/P  , where  f^fp)  «u(p>^^P  c^/p) . 

Then 

(2.11)  Cl(p,«)  S C(p,«)  i Cl((p,«)  . 

Proof;  See  Theorem  3.1  of  Werschulz  [76].  | 

Thus  we  have  bounds  on  the  complexity  of  using  pp  to  compute  an  approximate 
solution  satisfying  (2.4).  We  now  wish  to  consider  the  problem  of  optimality.  Define 

(2.12)  C*<«)  inf  {C(p,«);  pp  < ♦}  • 

We  are  interested  in  bounds  for  C*(e)  under  reasonable  assumptions  about  fj^  and  fy. 

We  first  suppose  that 

(2.13)  f|_(p)>0  and  fi/p)  > 0 if  p > 0 
and 

(2.14)  limptooft^P)  - ••"’ptoo  " • 

Assumption  (2.13)  is  that  there  is  no  method  whoso  cost  per  step  is  xero,  while  (2.14) 
essentially  means  that  the  "better"  a method  is  (i.e.,  the  higher  Its  order  is),  the  more 
we  should  expoct  to  pay  (or  its  use. 

Using  the  techniques  of  elementary  calculus,  we  find  that  a necessary  condition 


1 

t 


for  p to  minimize  C(_(  * ,e)  is  that 

(2.15)  • - Gl(p)  p^  fL'(p)  / (l(p)  I 
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similarly,  Cy( ' ,«)  takes  its  minimum  at  p pnly  if 

(2.16)  a - Gu<p)  p^  ^tp)  / fi/p)  . 

Sufficient  conditions  for  the  existence  and  uniqueness  of  solutions  to  \2.15)  and  (2.16) 
(i.e.,  for  well-defined  functional  inverses  of  (j|_  and  Gy)  which  actually  minimize  Cy(  * 
and  C|j(  ‘ ,«)  are  given  in 

Lemma  2.1;  Let  f|_  and  fy  be  as  above,  and  suppose  that 

(2.17)  Gl'(p)>0  if  Gl(p)>0  and  Gy'(p)  > 0 if  Gy(p)  > 0 . 

Then  G|_  and  Gy  have  respective  functional  inverses  P(_*  , py*  : R '*'**♦  F**  such  that 
for  all  p < R’*"*' 

(2.18)  CL*(a)  CL(P|.*(«),e)  S C^fp,®) 

and 

(2.19)  Cy*(a)  Cy(pu*(«),«)  S Cy(p,«) 

with  equality  in  (2.18)  or  (2.19)  if  and  only  if  p * Pl*(«)  or  p ■ py*(«),  respectively. 
Proof;  See  Theorem  2.1  and  Lemma  3.1  of  Werschulz  [76].  | 

We  call  Pl*(«)  (respectively,  Py*(o))  the  lower  (upper)  optimal  order.  Cl*(«) 
(respectively,  Cy*(«»  the  lower  (upper)  optimal  complexity,  and 

(2.20)  hL*<«)  hL(PL*<«),«)  (respectively,  hy*(«)  ;-  hy(py*(«),«)) 

the  lower  (upper)  optimal  step-size.  Combining  (2.1 1),  (2.12),  and  Lemma  2.1,  we  have 
Theorem  2.2; 

Cl*(o)  S C*(e)  S Cy*(«).  I 

We  next  describe  the  behavior  of  these  quantities  as  « increases  and  tends  to 
infinity. 

Theorem  2.3;  Let  fL  and  fy  be  as  in  Lemma  2.1.  Then  Pl*(«),  Py*(«),  Cl*(«),  and 
Cy*(«)  all  increase  monotonically  and  tend  to  infinity  with  a. 

Proof;  See  Theorems  2.2  and  3.3  of  Werschulz  [76].  | 
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Finally,  we  need  a restriction  of  the  problem  class  jDxV  to  "sufficiently  difficult* 
problems;  this  will  allow  us  to  determine  and  thus  establish  lower  bounds.  We  will 
assume  that 

(2.23)  i (Mi.h)P  if  h<  1 and  p< 

for  some  M(_  > 0 Independent  of  h and  p.  In  the  methods  we  study,  (2.23)  holds 
provided  all  sharp  upper  bounds  are  attained. 


||yil^jWi|lUilK^|iWWPWPP»,yP«>..-  ;|W.J.iirTIWiP»IWIWlill>i»Ml  mwi|.  ipMRMUPW  'iiJiu.jipwiiBippwwwii 
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3.  Taylor  Series  Methods 

The  class  of  Tavtor  series  methods  is  defined  by  expanding  x in  a truncated 
Taylor  series.  Thus  the  increment  function  is  given  by 

(3.1)  ^p(Xj,h)  v<*^>(Xj)  h*'  / (k+D! , 

where 

(3.2)  v<*'>(Xj)  (d/dt)*'  [v(x(t))l  | ^(t)  - Xj  • 

The  usual  method  of  computing  (3.2),  as  described  in  "classical"  numerical  analysis 
texts  such  as  Henrici  [62],  invokes  the  chain  rule.  This  quickly  leads  to  expressions  of 
horrifying  complexity;  for  this  reason,  most  texts  quickly  abandon  the  discussion  of 
high-order  Taylor  series  methods. 

We  are  interested  in  faster  algorithms  for  computing  ^p.  First,  we  address  the- 
problem  of  a lower  bound  for  the  combinatory  cost  d(p). 

Proposition  3.1;  There  exists  a constant  aj_  > 0 such  that  any  sequence  of 
algorithms  for  computing  must  satisfy 

(3.3)  d(p)  2 aLp'^  . 

Proof;  Any  algorithm  for  computing  fp  requires  the  information 
«p(v)  {D^v:  0s|(!|sp-  1)  . 

(We  use  the  standard  multi-index  notation  found  in  Friedman  [69].)  It  is  then  easy  to 
see  that  the  above  set  has  O(p^)  (as  p T oo)  distinct  elements,  which  are  (generally) 
Independent;  this  is  an  immediate  consequence  of  Problem  1 1 in  Chapter  I of  P6lya  and 
SzegO  [25].  Thus  (3.3)  gives  a linear  lower  bound.  | 

Note  that  the  constant  a|_  in  (3.3)  depends  on  N.  Since  we  are  treatir^  the  case 
where  N is  fixed  and  p is  allowed  to  vary,  we  will  not  indicate  this  dependence 
explicitly.  We  now  see  how  close  we  can  get  to  an  optimum  value  for  d(p). 


Theorem  3.1;  There  exists  a constant  ajj  > 0 such  that  the  combinatory  coat 
d(p)  of  computing  satisfies  the  bound 

(3.4)  d(p)  S ay  In  (p+e)  . 

Proof;  We  first  consider  the  case  N - 1.  Note  that  x(h)  is  the  zero  of 

(3.5)  F(z)  f J dt  / v(t)  - h . 

As  in  Brent  and  Rung  [76],  we  consider  the  formal  power  series 

P(s)  F(xq+s)  - F(xo> , 

where  s is  an  indeterminate.  Let  V be  the  power  series  reversion  of  P.  Adopting  the 
notation  of  Brent  and  Rung  [76],  we  see  that 

x(s)  - xo  + V(s)  - XQ  + Vp(s)  + CKsP^l)  . 

By  the  uniqueness  of  the  Taylor  coefficients  of  an  analytic  function,  we  see  that 

^p(xQ.h)  - h-^Vp(h). 

Since  the  number  Vp(h)  can  be  computed  in  (Mp  In  p)  operations  from  the  Taylor 
coefficients  of  v (by  Theorem  6.2  of  Brent  and  Rung  [76]),  the  result  for  N ■ 1 follows. 

For  N 2 2,  we  use  Newton’s  method  (Rail  [69])  applied  to  the  formal  power 
series  operator  P given  by 

(PyMs)  :»  y(s)  - xq  - Jq  v(y(r))  dr  } 

clearly,  the  formal  power  series  x(s)  is  the  zero  of  P.  The  algorithm  itself  is  defined 
recursively.  Let  a formal  power  series  x^p^(s)  satisfying 

X(p)(s)  - x(s)  + 0(sP*^ 

be  given.  Precompute 

(3.6)  w(s)  Jq  v(x(p)(r))  dr  - xq  - *(p)(s)  ♦ (Xs^P*^) , 

(3.7)  Q(s)  Vv(X(p)(s))  + 0(s2P*2) , 

and  let  u^q)(s)  0.  Then  set 


’‘(2p+l)<*>  *(p)<s>  ♦ “(prl)<*)  • 


li 


where 

(3.8)  U(^^j)(s)  !-  0(t)  U(|^)(r)  df  + w(s)  ♦ 0(s^^*^) , 0 S K S p . 

Following  the  proof  given  in  Rail  [69],  we  find  that 

X(2p+i)(»)  - x($)  ♦ (Xs2p*2)  . 

We  need  only  consider  the  cost  T(p,N)  of  computing  the  aerie*  x^p)(a)  in 
determining  d(p),  since  x(h)  may  be  recovered  from  the  formal  power  series  In  0(p) 
operations.  Clearly,  we  have  the  recursion 

(3.9)  T(2p+1,N)  S T(p,N)  + Tg  + T7  ♦ Tg  , 

where  is  the  cost  of  step  (3.m)  for  m ■ 6,  7,  8.  Let  COMP(p,N)  be  the  time  required 

to  find  the  first  p terms  of  the  formal  power  series  f(y^(s) y^ts)),  where  f,  yj, ...  , 

y|yj  are  formal  power  series,  and  yj, ... , y;^j  have  zero  constant  term.  Theorem  7.1  of 
Brent  and  Kung  [76]  states  that 

C0MP(p,2)  - (Xp2|np), 
and  it  is  easy  to  show  that  for  any  N < 1**, 

COMP(p,N+l)  - CXp  COMP(p,N))  . 

Thus  for  N 2 2,  we  have 

(3.10)  C0MP(p,N)  - 0(0*^10  0), 
and  so  we  see  that 

Tg  + T7  - 0((2p+l)N|np). 

Finally,  let  MULT(p)  be  as  in  Brent  and  Kung  [76J  we  see  that 

Tg  - (p+1)  [n2  MULT(2p+l)  ♦0(p)]  - 0((2p+l)2  In  p) 
if  Fast  Fourier  Transform  multiplication  (Borodin  and  Munro  [75])  is  used.  Since  N 2 2, 
we  have 

(3.1 1)  Tg  ♦ T7  ♦ Tg  - CX(2p+l)N  In  p) , 


and  so  (3.9)  and  (3.11)  imply  that 
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T(p,N>  - CXp*^  tn  p) , 

which  completes  the  proof.  | 

(Note  that  the  second  algorithm  is  Inferior  to  the  first  algorithm  when  applied  to 
the  scalar  case  N - 1,  where  we  find  that  the  second  algorithm  requires  (Xp^  In  p) 
arithmetic  operations.) 

We  now  determine  bounds  on  C(p,«).  First,  consider  lower  bounds.  Clearl/, 
there  exists  e|_(v)  2 0 such  that 

(3.12)  e(D^Vj)  i sl(v)  (1  S i S n,  |ff|  € I*)  . 

Since  9tp(v)  has  0(p^)  elements,  there  exists  a constant  b|_  > 0 such  that 

(3.13)  e(Wp(v))  i bL  Ol(v)  pN  . 

From  (3.3)  and  (3.13),  we  have  a lower -bound  cost  per  step  of 

(3.14)  Cl(p)  - [sl  + bL  6l(v)]  pN  . 

This  leads  to 

Theorem  Cl(p,«)  ■ Ml  [aL  + bL  eL(v)J  p*'*  e*/P  . 

Proof;  This  is  an  immediate  consequence  of  (2.23)  and  (3.14).  | 

Note  that  (l(p)  :>  MlCl(p)  satisfies  the  conditions  of  Lemma  2.1.  Thus,  the 
optimality  theory  of  Section  2 holds.  In  particular,  we  have 
Theorem  Cl*(«)  - Ml  [sl  + bL  ol(v)]  (e/N)^  . 

Proof;  From  (2.18)  and  (3.14),  we  find  that  Gl(p)  ■ Np,  so  that 
PL*(o)  - «/N  and  hL*(«)  - (MloN)-!  . 

The  result  follows  by  letting  p • Pl*(«)  in  the  definition  of  Cl(p,«).  | 

However,  recall  that  we  assumed  that  the  non-identical  mixed  partial  derivatives 
of  V are  independent.  There  are  a number  of  systems  for  which  this  is  not  true  (for 
instance,  constant  coefficient  linear  systems);  for  such  systems,  >t  is  clear  that  we  may 
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be  able  to  use  the  extra  information  of  non-independence  to  find  algorithms  that  are 
faster  than  the  lower  bounds  given  above.  However,  we  will  ignore  this  case  and  only 
consider  the  problem  for  a "general"  function  v. 

Next,  we  turn  to  upper  bounds  on  the  complexity.  Theorem  3.1  tells  us  how  to 
combine  the  necessary  information  to  get  the  solution  at  a new  grid-point}  we  need 
only  measure  the  cost  of  getting  the  information.  So,  let 

e<'^)(v)  - max  {e(D^vj):  1 S i s |fl|  - k)  . 

Using  the  result  in  P6lya  and  SzegS  [25],  we  see  that 

(3.15)  e(9lp(v»  S N e<^v)  (N*K-D!  / [k!(N-D!]  . 

Unfortunately,  the  right-hand  side  of  (3.15)  does  not  fit  our  general  model,  so  we  must 
assume  that  we  know  how  e^*^\v)  changes  as  k increases.  Wg,  wjll  consider  the  case 
where  the  cost  of  derivative  evaluation  is  bounded;  that  is.  we  will  assume  that 

(3.16)  e^*'\v)  ^ eu(v) 

ffir  some  eu(v)  independent  gf  k.  Other  cases  (e.g.,  e^*^^(v)  - Otk"*)  for  some  m > 0) 
may  be  analyzed  in  a similar  manner;  of  course,  they  will  give  different  results.  By 
(3.15)  and  (3.16),  there  is  a by  > 0 such  that 

(3.17)  e(92p(v))  S byeu(v)p*'*. 

From  (3.4)  and  (3.17),  we  have  an  upper-bound  cost  per  step  of 

(3.18)  Cl)(p)  - ay  p^  In  (p4e)  ♦ by  ey(v)p*'* . 

This  leads  to 

Theorem  3.4:  There  exists  an  My  > 0 such  that 

Cy(p,«)  - My  [ay  p*'*  In  (p+e)  ♦ by  ey(v)p^]  e*/P  . 

Proof;  By  Cauchy's  Integral  Theorem  (AhIfors  [66],  pg.  122),  there  exists  a 
B > 0 such  that 

|||x<''*i>|||/(k+l)l  i B^  , 
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where  we  define 

(3.19)  lllylil  J-  "*«  t < I 

for  any  y;  1 -»  R^.  Thus  by  Section  3.3-3  of  Henrici  [621  we  see  that  a UpaeNt* 
constant  for  pp  in  # j is  given  by 

|||x<'‘^l>|||  h'^  / (K+D!  S SflJ(Bh)''  S L (1  - Bho)'l  , 
provided  that  h S hQ  < By  Section  3.3-2  and  3.3-4  of  Henrici  [621  there  exiata  an 
My  > 0 such  that 

^ (K^h)P. 

The  result  now  follows  from  Theorem  4.1  and  (3.18).  | 

We  are  now  ready  to  consider  the  optimal  p for  C(j(p,«). 

Theorem  3.5; 

(1.)  For  all  a > 0,  there  exists  py*(«)  such  that  (2.19)  holds. 

(2.)  Py*(«)  increases  monotonically  with  a,  and 

py*(«)  a/N  as  a t 00  . 

(3.)  Cy*(a)  increases  monotonically  with  a,  and 

Cy*(a)  ~ My  ay  (e/N)*'*  a^  In  a as  a t oo  . 

(4.)  hy*(a)  ~ (My  eN)“l  as  a T oo  . 

Proof;  Clearly  Cy  satisfies  (2.13)  and  (2.14).  Now  write 

Gy(p)  • Gi(p)  ♦ G2(p) » 

where 

Gj(p)  - Np  and  G2(p)  - rp^/02(p)l 

here  we  set 

02(p)  !■  (p+e)  [(p+e)  In  (p+e)  ♦ 1]  and  r ;-  ay/[byey(v)]  . 

We  see  immediately  that  Gj  satisfies  (2.17);  a straightforward  calculation  shows  that 
G2'(p)  - 9 [0(p)]'2  {xp  [In  (p4e)]  - 1]  ♦ 2e[a  In  (p4e)  ♦ 1]} , 
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so  that  G2'(p)  > 0 for  p > 0.  Thus  63  satisfies  (2.17).  which  shows  that  Gy  satisfies 
(2.17).  Hence  py*  and  Cy*  behave  as  described  in  Theorem  2.2. 

Since  py*(«)  goes  to  infinity  with  e,  we  see  that 

« ■ Gi/Pu*^*^  N Pu*(«^  ♦ PU*<*^  / ^ PU*<*^  * 

which  gives  the  asymptotic  estimate  in  (2.).  The  rest  of  the  The  rem  follows  from  this 
estimate.  | 

Unfortunately,  the  estimates  given  above  are  only  asymptotic  as  • t ooi  this  will 
be  typical,  since  many  of  the  equations  to  be  solved  involve  products  of  logarithmic 
and  polynomial  terms,  and  thus  cannot  be  solved  exactly.  On  the  other  hand,  these 
asymptotic  expressions  are  sufficient  for  our  purposes,  since  they  describe  how 
quickly  py*(«)  and  Cy*(«)  increase  with  a. 

Note  that  as  a tends  to  infinity,  Cy*(«)  becomes  independent  of  ey(v),  which 
measures  how  hard  it  is  to  evaluate  the  derivatives  of  v;  this  is  because  the 
combinatory  cost  eventually  overwhelms  the  informational  cost.  This  Kind  of  behavior 
will  be  typical  of  the  complexity  analyses  in  this  paper.  Finally,  note  that  the  bound 

(3.20)  CL*(a)  - 0(eN)  s C%)  i 0(«'^lno)  - Cy*(e)  as  « t 00 
implies  that 

Cy%)  / Cl*(«)  - 0(ln  e)  as  « t 00  j 

this  indicates  the  gap  in  our  knowledge  of  the  complexity  of  solving  (1.1)  via  Taylor 
series  methods. 


‘1 

i'] 


I 


( 
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4.  Linear  Runge-Kutta  Methods 


For  many  functions  v,  caculation  of  tha  derivativos  requirad  by  Taylor  seriaa 
methods  is  prohibitivaly  axpansiva.  For  this  raason,  wa  ara  intarastad  in  mathods 
which  use  information  that  is  somewhat  more  raadily  availabla.  In  particular,  wa  will 
considar  mathods  that  usa  only  avaluations  of  v,  combinad  in  a higMy  atructurad 
manner.  Wa  say  that  is  a class  of  linear  Runae-Kutta  methods  (abbraviatad,  IRK 
methods”)  if  each  increment  function  m>y  be  written  in  tha  form 


:■  Xji  K| 


where 


(4.2)  k,  v(Xj  ♦ h kj)  for  0 i I S s - 1 , 

the  integer  s - s(p)  is  said  to  be  the  number  of  stages  of  pp;  the  number  of  stages  is 
equal  to  the  number  of  times  tha  vector  furKtion  v must  be  evaluated.  (In  order  to 
simplify  notation,  we  will  not  explicitly  indicate  the  dependence  of  Af j and  kj  on  p.)  The 
method  ep  defined  by  (4.1)  and  (4.2)  is  explicit  in  that  k|  depends  only  on  k^  , k|.|{ 
see  Butcher  [64]  for  a discussion  of  semi-exolicit  and  implicit  mathods.  (Wa  use  the 
adjective  ”(inear”  to  distinguish  these  methods  from  "nonlinear  Runge-Kutta  methods," 
which  were  first  proposed  in  Brent  [74].) 

Since  the  function  ep  is  (in  practice)  always  evaluated  by  using  the  obvious 
algorithm  suggested  by  its  definition,  we  shall  identify  an  algorithm  for  evaluating 
with  ep  itself.  Thus  the  problem  of  finding  the  best  algorithm  for  evaluating  ^p  in 
^LRK  equivalent  to  the  problem  of  finding  the  best  basic  sequence  of  LRK  methods 
possible.  This  is  related  to  the  problem  of  finding  the  smallest  value  of  s(p)  such  that 


ep  has  order  p.  This  minimal  value  is  given  by 


unknown 


For  methods  of  order  greeter  then  seven,  e gep  develops.  For  instence,  eighth-order 
methods  with  eleven  steges  exist,  end  it  is  known  thet  eny  eighth-order  method 
requires  et  leest  ten  steges.  For  erbitrery  p 2 8,  the  best  bounds  known  for  the 
optimum  velue  of  s(p)  ere 


(4.4)  p ♦ #(p)  s s(p)  i ip‘  - 7p  ♦ 14)  / 2 , 

where  #<p)  2 c In  p for  ell  sufficiently  lerge  p (for  some  c > 0).  The  lower  bound  is 
given  in  Butcher  [75]i  the  proof  is  quite  involved,  end  the  result  is  not  much  better 
then  the  Triviel"  lower  bound  s(p)  2 p (Hkidmersh  [74X  pege  84).  A dess  #cvRK 
methods  such  thet  pp  requires  only  (p^  - 7p  ♦ 14)  / 2 steges  is  given  in  Cooper  end 
Verner  [72J 

tMe  first  consider  lower  bounds  on  the  complexity  C(ppi)  using  LRK  methods. 
The  "trivisT  lower  bound  s(p)  * ^0  will  be  used,  since  the  term  #(p)  will  be  smell  when 
p is  smell  end  wHI  not  effect  the  esymptotic  behavior  of  optimel  order  end  complexity 
for  p lerge.  It  is  known  (Butcher  [84])  thet  et  leest  0(p^)  of  the  subdiegonei  elements 
of  the  metrix  A (whose  elements  ere  the  k||  in  (4.2))  must  be  non-xero  in  order  for  A 
to  define  e p^-order  method.  Thus  there  exists  Sf^  > 0 such  that 

(4.5)  d(p)  2 a|^  t 
since  s(p)  2 p,  we  see  thet 


where  we  now  write 


? 

I 


Thu*  (4.6)  and  (4.6)  show  that  a lowar  bound  on  tha  coat  par  atop  for  h givan  by 
(4.7)  C|_(p)  - *L  p2  + N 0|_(v)  p . 


IbSfiCftmM: 


Cl(P^)  - ML[aLp2  + N*L(v)p]a«/P  . 


Proof!  Thia  follows  immadiataly  from  (2.23)  and  (4.7).  | 

It  is  claar  that  fj^tp)  :•  Ml  f*L  ♦ N ol(v)  p]  a*^P  satisflas  (2.13)  and  (2.14). 
Wa  claim  that  fj^  yields  a G(_  satisfying  (2.17).  Indeed,  write 

ft<P^  ■ U^P^^2^P^' 


where 


fj(p)  Ml  *L  P 


f2(p)  P * • where  r :•  N e^tv)  / sl  . 

Clearly  fj  yields  a Gj  satisfying  (2.17).  Since  (2  is  a linear  polynomial  with  a negative 
zero,  it  may  be  shown  that  (2  yields  a G2  satisfying  (2.17).  Thus  f|_  yields  a G|_ 
satisfying  (2.17h  in  fact,  we  have 

(4.8)  G^tp)  - Gi(p)*G2(p)  - p[l  ♦(!  • 

This  leads  us  to 

Theorem  4.2: 

Cl*(«)  [Ml  sl  e^  / 4]  ar^  as  « T os  . 

Proof:  From  (4.8),  we  see  that  Gl(p)  2 p as  p t oo.  Since  (2.13),  (2.14),  and 
(2.17)  hold,  Pl*(«)  tends  to  infinity  with  «.  Thus 

• - Gl<Pl*<«))  2 Pl*(«)  as  « t 00  , 

i.e.,  PL*(e)  e/2  as  « t 00.  The  result  now  follows  from  Theorem  4.1.  | 

Wa  now  turn  to  upper  bounds  on  complexity.  The  class  #cvrk  derived  in 
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Cooper  and  Verner  [72]  has  two  deficiencies,  the  first  of  which  is  that  rw  uniform 
upper  bound  on  e|_|j(^p,h)  is  Known  for  #CVRK’  addition,  the  combinatory  cost  for 
this  class  of  methods  is  CKp^)  as  p t oa.  Instead,  we  turn  to  the  basic  sequence  #CRK 
discussed  in  the  Appendix.  There,  we  prove  that  there  is  an  > 0 such  that 

(4.9)  ^ <P  ♦ ®>  » 

provided  h S hp,  where  hp  - 0«ln  p)“^)  as  p t oo.  Furthermore,  there  are  a large 
number  of  extra  zeros  in  the  matrix  A for  pp  < Using  the  notation  of  the 

Appendix,  we  see  that  the  number  of  non-zero  entries  in  A is 

(I  ■ Cl' * " 

- p3/3  - p2/2  ♦ 7p/6 

S p3/3  + 2p2/3 

for  p € Z Finally,  note  that  the  number  of  stages  s(p)  required  for  pp  ( #^RK 

(4.10)  s(p)  - l(p2  - 2p  ♦ 4)/2J  i p^l2  * p 


for  p < 1 which  shows  that  the  number  of  stages  required  for  a p^^-order  method 
in  asymptotically  equals  the  number  requires  for  a p^^-order  method  in  #CVRK' 

Thus  (considering  the  combinatory  costs),  the  class  ^qvpk  actually  costs  more  per 
step  than  does  #CRK'  >S*^oring  the  combinatory  costs  would  have  caused  us  to  reach 
the  opposite  conclusion. 

First,  we  look  at  the  cost  per  step.  By  (4.10), we  see  that 

(4. 1 1 ) e( Wp(v))  S - (p2  ♦ p)  N ei/v) , 

where 

eu(v)  max  e(vj)  . 

Since  we  are  using  it  is  easy  to  see  that  there  is  a b(j  IC  2/3  such  that 

(4.12)  tKp)  i (p3/3  ♦ bu  P^>  • 2N  . 


i 

! 


Combining  (4.11)  and  (4.12),  we  see  that  the  total  combinatory  cost  per  step  it 
bounded  by 

(4.13)  cu(p)  - N [ 2p3/3  ♦ p2  ♦ ^2  P 1 • 

where 

fi\  5“  ®i/'')  / 2 ♦ 2 bjj  and  $2  •"  ^ • 

Using  (4.9)  and  (4.13)  gives 
Theorem  4.3; 

Cu(p,«)  - Mu  N [2p3/3  ♦ ill  p2  ♦ II2  P 5 <P  * • ■ 

Now  we  looK  at  the  optimality  theory  for  the  upper  bound. 

Theorem  4.4.; 

(1.)  For  ali  « > 0,  there  exists  pu*(«)  such  that  (2.19)  holds. 

(2.)  Pu*(«)  increases  monotonically  with  % and 
PU*(«)  «/3  as  « t 00  . 

(3.)  Cu*(«)  increases  monotonically  with  «,  and 

Cu*(«)  ^ [ 2 Mu  N e^  / 81  ] In  e as  e t 00. 

(4.)  hu*(«)  ^ (Mue^lne)"^  as  • t 00  . 

Proof:  We  write 

fu(p)  Mu  In  (p  ♦ e)  cu(p) 

in  the  form 

fu(p)  • fi<p)f2(p)» 


^l(p)  ■ MuNpln(p  + e)  and  f2<p)  ■ 2pV3  ♦ p ♦ ^2  • 

It  is  clear  that  f 1 satisfies  the  hypotheses  of  Lemma  2.1.  Now  we  consider  (2-  Clearly 
fo  has  no  positive  zeros;  it  may  be  seen  that  the  condition  bu  i 2/3  implies  that  (2  has 


a positive  discriminant  and  hence  has  no  complex  r ^ots.  Thus  f2  has  only  negative 
rootsi  one  may  then  show  that  this  guarantees  that  f2  satisfies  the  hypotheses  of 
Lemma  2.1.  Thus,  the  same  may  be  said  for  f * fj  f2  • 

Thus  P|j*  and  Cj*  behave  as  described  in  (1.)  of  Theorem  2.3.  We  also  see  that 
G(j(p)  « 3 p as  p t 00.  Thus  the  estimate  in  (2.)  holds,  from  which  we  get  the  estimates 
in  (3.)  and  (4.).  | 

So  in  the  class  of  linear  Runge-Kutta  methods,  we  find  that 
(4.14)  CL*(e>  - 0(e2)  s C*(e)  i - (Xe^ 

as  « tends  to  infinity:  hence,  the  ratio 

Cu•(e)/CL^e)  - CKe  In  e) 

indicates  the  gap  in  our  Knowledge  of  the  complexity  of  linear  Runge-Kutta  methods. 


5.  Comparison  of  the  Methods 


We  now  wish  to  compare  the  classes  of  Taylor  series  methods  and  LRK  methods. 
Write  » ^T*  • ^T*  (•’•*P«c^'valy,  Cjj^lRK*  » ^LRK*^  ^U* 

and  C*  in  the  class  fj  (respectively,  the  class  #i.pk)*  Since  we  have  only 
asymptotic  expressions  for  these  quantities,  we  are  forced  to  use  an  asymptotic 
comparison.  If  f,  g ; F**  R'*’*  satisfy  lim  f(«)  ■ lim  g(«)  ■ ♦<»,  we  wHI 

write 

(5.1)  f < g iff  f(«)  ■ o(g(«))  as  « t 00  { 

we  say  that  f is  asymptotically  less  than  g.  If  f < g,  there  is  an  > 0 such  that 
f(«)  < g(a)  for  a > oq,  so  there  is  a non-asymptotic  interpretation  of  the  order 
relation  <.  Thus  if  f and  g are  cost  functions,  the  statement  "f  < g"  implies  that  the 
method  whose  cost  is  given  by  f is  "better"  (i.e.,  cheaper)  than  the  method  whose  cost 
is  given  by  g,  for  t sufficiently  small.  Using  the  results  of  (3.20)  and  (4.14),  we  then 
have  the  following 

Theorem  5.1;  Suppose  that  (3.16)  holds. 


(1.) 

If  N - 1,  then 

Cu,T*  ^ • 

(2.) 

If  N - 2,  then 

^U,T*  ^ ^LRK*  • 

(3.) 

If  N - 3,  then 

Cu^T*<«>  • 

and 

Cuj.RK*<«>  • 

as  • t 00. 

(4.)  If  N ^ 4,  then  C^lrk*  < • ■ 


be  true.  As  sn  immediate  coroliary  to  the  above  theorem,  we  have 
Theorem  5.2; 

(1.)  If  N - 1 and  (3.16)  holds,  then  Ci-*  < C|,pK*  . 

(2.)  If  N i 4,  then  C^rk*  < Cj*  . I 

So  if  the  derivatives  of  v are  cheap  to  evaluate,  we  see  that  the  best  Taylor 
series  method  known  is  better  than  the  best  linear  Runge-Kutta  method  possible  for 
the  scalar  case  N > 1;  but  if  N ^ 4,  the  best  linear  Runge-Kutta  method  known  is  better 
than  the  best  Taylor  series  method  possible. 


24 


Appendix:  Error  Bounds  for  a Sequence  of  LRK  Methods 


I ! 


In  this  Appendix,  we  describe  a subclass  of  a class  of  linear  Runge-Kutta  ("LRK“) 
methods  due  to  Cooper  [69].  We  shall  first  prove  the  following 

Theorem  A.!;  There  is  a basic  sequence  of  LRK  methods  such  that 

(1.)  Each  iPp  i *cpK' 

s(p)  :•  (p2  - p ♦ 2)  / 2 
evaluations  of  v per  step. 

(2.)  There  exists  an  My  > 0 such  that 

(A.1)  In  (p+e)  h)P 

for  h S hp  - CXdn  p)*^). 

We  use  the  notation  of  Cooper  and  Verner  [72].  Let  p < Z'*"*’  bo  given;  define 


p:  Z n [0,  p]  -»  Z * by 


(A.2)  p(j> 


^h-0  ^ ^ 'M  P 

s if  j - P • 


where  we  write  "s"  for  "s(p)"  as  defined  above.  Next,  a set  {{q  , ... , of  integers  is 
defined  by  picking  (q  :■  p,  and  setting  (j  (i  0)  to  be  the  unique  integer  in  [1,  p] 


satisfying 


p(fe  - 1)  < i S p(fc)  . 


We  now  pick  Uq  , ... , u^  < I satisfying 


Uq  • 0 , Ug  - 1 , Uj  0 if  i 0 


({j  " (j  and  i i*  j)  implies  Uj  ^ Uj 


Finally,  we  pick  a matrix  of  coefficients  A :■  {kjj:  0 S j S i-l,  1 S i S s)  such  that 
(A.6)  Ajj  - 0 if  Ij  < Ij  - 1 (1  S i,  j S s) 
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and 

(A.7)  Xjj  U|»  - (0  S r S tj  - 1, 1 S I S a)  . 

Cooper  and  Verner  [72]  point  out  that  these  conditions  may  always  be  fulfilledi  the 
resulting  A defines  a p^^-order  LRK  method  with  s stages. 

We  are  interested  in  a choice  of  uq u,  which  will  give  a small  error 

coefficient.  To  this  end,  we  wiil  choose 

(A.8)  {Uj5  - {(l+X|^„)/2:  ISKSn)  (1  S n S p - 1) , 

where  , ...  , are  the  zeros  of  the  Jacobi  polynomial 
SzegS  [59]).  Since  ihese  zeros  are  distinct  and  lie  in  [-1,  1] , conditions  (A.A)  and  (A.5) 
may  be  satisfied. 

Now  we  are  abie  to  exhibit  a solution  to  the  i^^  system  In  (A.7).  First,  note  that 
the  equation  for  r - 0 may  be  separated  from  the  others,  since  uq  - 0.  Setting 

n - I , 

we  see  that 

(A.9)  Xjo  - Uj  - - Uj  - I { Xjj:  j < i and  a n } , 

the  last  by  (A.6).  We  wish  to  determine  the  nonzero  Xjj,  i.e.,  those  Xjj  for  which  |j  i n 
and  j < i.  So  setting 

Xjj  - 0 unless  j < {ji  , ... , ]„} . 
we  see  that  the  remaining  Xjj  are  the  solution  of  the  system 

Thus  the  Xjj^  gra  [he  weights  [or  gn  interpolatory  quadrature  formula  fin  [0,  U|]  jjtiHi 

abscissae  u:  u:  . From  the  usual  exoression  for  such  weights  and  (A.6),  we  see 

>l  In 

that 

“ Pikn  [2Pn'(cos  d^„)]‘l  ^ ^ ‘ . 

where  X|^^  - cos  (1  S K $ n). 
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Lftmmi  AJ.S  HjKn  “ n)  «8  n t 00. 

Proof;  Since  the  zeros  of  ere  symmetric  about  the  origin,  we  may  eeeume 
that  0 < #|^n  S w/2.  Using  (8.9.2)  of  Szeg6  [591 

Pikn  ■ 0(K5/2n-3)  [P„(cos  #)  / (cos  # - cos  #Kn>l  • 

Cltt  1:  ♦l.n+1  ^ *i,n+l  ^ ^K,n*l/2-  ^ 

[#l„/2.  #j  n+j]  , since  Theorem  15.4  of  Szegft  [59]  proves  that 

0(k5/2n-3)  [I  J'  I ♦ I J II  - 0(n-l)  . 

(Here  the  integrand  is  the  same  as  in  the  preceding  integral.)  But  the  proof  of  (15.4.12) 


in  SzegS  [59]  extends  almost  immediately  to  a proof  that  the  remaining  integral  is 
(Xk"2n),  since  (15.4.12)  is  proved  by  order-of-magnitude  estimates.  Thus  «■{)„  • 
CXn"^)  - 0(n“^  In  n)  for  Case  1. 

Casft  2:  ^k.n*l/2  ^ ^i,n*l  s 3tfk,n*l/2-  Vife  consider  the  integral  over 
[d^„/2,  ^j,n+l]  » Szegfi  [59]  shows  that 

0(k5/2n-3)  If'  I - 0(n-l)  . 

•^'kn/2 

As  in  (15.4.13)  of  Szegft  [591  we  have 

♦ >2  • 


Here 


with 


«kn/2 

II  5-  0(d)  sin  d dd  , 

^ ^^kn/2 


D(d)  [cos  (Nd  ♦ y)  - cos  (Nd,^„  ♦ t)]  / [cos  d - cos  d|^„] , 


where  N ;■  n + 3/2  and  y :■  -3w/4,  and 

Ig  J R„(d,d^„)  sin  d dd  - 0(nk-3/2) , 
with  the  remainder  term  in  (8.8.2)  of  Szegfi  [591  Unfortunately,  the  proof  that 
(15.4.14)  of  SzegO  [59]  is  bounded  does  Qfit  extend  to  a proof  that  Ij  is  bounded, 
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sine©  the  proof  of  th©  form©r  r©qulr©s  that  th©  int©rv©l  of  intogration  b©  aymm©trlc 
about  How©v©r,  it  it  straightforward  to  vorify  that 

>1  ■ 0(1)  JJ^^|sinN#/#|dd  - O(lnn)  . 

Thus  ■ (Xn'^K  In  n)  - 0(n"^  In  n)  for  Cat©  2. 

Cas©  3i  3#K,n4.1  i ^ 3*/^-  consld©r  th©  int©gral  ov©r 

I Sz©gO  [59]  prov©s  that 

CKk6/V3)  I I . an-l)  . 

But  th©  proof  of  (15.4.19)  in  Szegfl  [59]  extonds  to  prov©  that  th©  r©maining  integral  is 
0(|^”5/2n)  (as  in  Cas©  1).  Thus  « 0(n"^)  • 0(n"^  In  n)  for  Cas©  3. 

Case  4:  3w/4  S S *n*l,n+l*  consider  th©  integral  over 

[3»/4,  #j^n+l3  • Sz©g5  [59]  shows  that 

0(k5/2n-3)  I I - (Xn-h  . 

As  in  Cases  1 and  3,  th©  proof  of  th©  above  may  b©  extended  to  prov©  a similar  bound 
on  th©  integral  of  interest.  Thus  ii^n  ■ 0(n~h  “ 0(n"^  In  n)  in  Cas©  4,  completing  the 
proof  of  th©  Lemma.  | 

Thus  (A.9)  and  Lemma  A.1  show  the  existence  of  a X > 0 such  that 
(A.10)  sj;{j  IXjjl  S X In  (tj  ♦ ©) ! 

her©  X is  independent  of  p.  Moreover,  the  result  for  the  cas©  i - s may  be  sharpened. 
We  see  that  X-:  2 0,  since  the  u:  for  the  s*^  system  in  (A.7)  are  the  abscissae  for 

I 


Lobatto  quadrature.  Thus 
(A.11) 


1^)1  ■ *-:i  m • >• 


the  consistency  condition  In  the  last  equality  being  a consequence  of  (A.7)  with  t ■ 0. 
Proof  2t  Theorem  A.l;  As  in  Cooper  and  Verner  [72],  wo  define 


ij  i(Ujh)  - Kj 
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for  0 S i S St  not*  that  Iq  " *0  “ **•  computed  ipproxlmetton  to  x(h)i 


(A.12) 


h-1  ||K<h)  - r(h)||  - II  h-1  [x(h)  - x(0)]  - X,,  Kj  || 

s in.li  ♦ P?:;5  


i l|S,||  ♦ max  . p_i  lltill , 

the  last  by  (A.6)  and  (A.1 1).  By  the  analyticity  of  x,  there  is  an  Aj  > 0 such  that 
fit  j-  h-1  II  x(Ujh)  - xJLo  (Ujh)V»>(0)  / Tt  II  s (Aj  h)*^ 
and 

r,i  II  «Ujh)  - (ujh)»  4<»>0)  / r!  II  S (Aj  h)^  , 
so  that  the  definition  of  Ij  gives 


(A.13) 


ll»ill  i ^i  ♦ *;■-!)  INil  tij 

i (Ai  h>^'  ♦ IXjjl  <Ai  h)^ 


S (A2hr 

for  a suitable  A£  > 0.  Thus  (A.12)  becomes 

(A.14)  h"^  |(x(h)  - 2(h)||  i (A2  h)P  ♦ max  , p.j  HsjH  . 

We  now  use  Lemma  1.1  of  Cooper  and  Verner  [72]  and  (A.6)  to  fifKl  that  if  L is  a 
Lipschitz  constant  for  v,  then  there  exists  A3  > 0 such  that 


Ilf, II  i hL  ||6,||  ♦ hL  IXjjl  max  j ||fj|| 

i (A3  h)^'^^  ♦ (A3  h)  In  ((,  ♦ e)  max  j ||ij|| , 

the  last  by  (A.10)  and  (A.13)t  here,  the  maximum  is  taken  over  all  J < I such  that 
Ij  2 - 1.  A straightforward  induction  shows  that  if  (1  ♦ In  2)  A3  h < 1,  then 

tisjll  S <A4  In  ({,  ♦ e)  h)*'^^ 
for  a suitable  A4  > 0.  Combining  this  with  (A.14),  we  find 

(A.15)  h"l  ||x(h)  - z(h)|l  i (A5  In  (p+e)  h)P  , 
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I 

t 


l«»K'  • - • 


the  desired  bound  for  the  local  error  for  a single  unit  step. 

To  extend  (A.  15)  to  a global  error  result,  we  must  look  at  the  Lipschitz  constants 
for  the  increment  functions.  Let  L be  a bound  on  ||t7v||,  and  write  "Vpptyih)"  to 
indicate  gradient  with  respect  to  the  vector  variable  y.  Now 
|IWp(y,h)||  s Ik^il  max  o^iss-l 

- Osiss-1  » 

where  we  write  1ij(y,h)"  to  indicate  the  dependence  of  kj  upon  y and  h.  By  the 
definition  of  kj(y,h),  we  find 

Vkj(y.h)  - 7v(u)  [l^xN  ♦ h xj;},  Xjj  7kj(y.h)] , 
where  u ;■  y + h x|]|^  Xjj  kj(y,h)  and  is  an  NxN  identity  matrix.  Taking  norms  In 
the  above  gives  the  result 

fj  S LX  ♦ hLX  [ In  (|j+e)  max  { f j < i and  |j  2 - 1 } ] • 

where  fj  :■  t|ykj<y,h)f|.  Writing  Xp  for  the  Lipschitz  constant  for  pp,  it  is  easy  to  see 
that  (A.  16)  and  the  above  inequality  imply 

Xp  S Xj*^^  (hLX)i  nj'fj  In  (p+e-k) , 

which  is  bounded  for  all  p,  provided  that  h s hp  < (LX  In  (p+e))~^.  Thus  (A.1)  follows 
from  this  result,  (A.  15),  and  Theorem  3.3  of  Henrici  [62].  | 

The  value  for  s(p)  indicated  in  Theorem  A.1  may  be  improved  somewhat  by 
noting  that  since  we  are  using  a Lobatto  quadrature,  higher  order  may  be  expected 
with  fewer  steps.  Indeed,  if  we  use  the  strategy  outlined  in  the  comments  following 
Theorem  4 of  Cooper  and  Verner  [72],  we  have 

Theorem  A.2;  There  exists  a basic  sequence  4cRK  methods  such  that 

(A.1)  holds  and  ^p  requires 

s(p)  :•  Kp*  - 2p  ♦ 4)  / 2 j 

evaluatiorts  of  v per  step.  | 
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